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We show that density functional theory within the RPA (random phase approximation for the 
exchange-correlation energy) provides a correct description of bond dissociation in H2 in a spin- 
restricted Kohn-Sham formalism, i.e. without artificial symmetry breaking. We present accurate 
adiabatic connection curves both at equilibrium and beyond the Coulson-Fisher point. The strong 
curvature at large bond length implies important static (left-right) correlation, justifying modern 
hybrid functional constructions but also demonstrating their limitations. Although exact at infinite 
and accurate around the equilibrium bond length, the RPA dissociation curve displays unphysical 
repulsion at larger but finite bond lengths. Going beyond the RPA by including the exact exchange 
kernel (RPA+X), we find a similar repulsion. We argue that this deficiency is due to the absence of 
double excitations in adiabatic linear response theory. Further analyzing the H2 dissociation limit 
we show that the RPA+X is not size-consistent, in contrast to the RPA. 

PACS numbers: 31.15.Ew, 31.25.Eb, 31.25.Nj 
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I. INTRODUCTION 

Density functional theory [1-3] (DFT) has proved to 
be a powerful method for calculating (and analyzing) 
the ground-state properties of molecular and condensed 
matter. In its standard Kohn-Sham (KS) form, the 
density n(r) and total energy are constructed from the 
self-consistent solution of one-electron equations where, 
in practice, the exchange-correlation (XC) energy func- 
tional i5 xc [n] must be approximated. Already the simple 
local-density approximation (LDA) and, more so, gener- 
alized gradient approximations (GGAs) can yield a re- 
markably realistic description of chemical bonds in solid 
and molecular systems. The state of the art is presently 
set by hybrid functionals [4, 5] that admix a fraction of 
the exact (Fock) exchange energy with GGA exchange. 
Achieving on the average nearly chemical accuracy for 
bond dissociation energies, they rival much more de- 
manding post-Hartree-Fock configuration interaction or 
coupled cluster methods. 

However there remain well-known conceptual limita- 
tions, with a clear practical significance, as exemplified 
by the following paradigm situations, (i) Long-ranged 
Coulomb correlations between non overlapping systems 
are not included in local functionals such as the LDA or 
GGAs (and thus also hybrids) but require fully nonlo- 
cal XC energy functionals [6-8] . Indeed LDA and GGAs 
perform at best erratically for van der Waals bonded sys- 
tems [9, 10]. (ii) Scaling to the high-density or weakly- 
interacting regime, hybrid functionals do not properly 
recover the exact KS exchange energy [11]. This fail- 
ure prominently concerns odd electron bonds, as ex- 
emplified by the molecule, where 100% exact ex- 
change mixing and zero correlation energy would be 



needed [12, 13]. (Hi) In systems with significant static 
(non-dynamical) correlation, LDA, GGA, and thus hy- 
brid functionals underestimate the magnitude of the cor- 
relation energy [14, 15]. This becomes particularly prob- 
lematic for the dissociation of electron pair bonds where 
near degeneracy effects arise in the molecular wavefunc- 
tion. A famous example is the dissociating H2 molecule: 
the proper (singlet) KS ground state at larger bond 
lengths has much too high total energy for such func- 
tionals. Usually one works around the problem: per- 
forming a spin-unrestricted calculation, a second solu- 
tion with reasonable, lower energy is obtained. This suc- 
ceeds because exchange functionals can mimic the static 
(i.e. long-ranged left-right) correlation and thus com- 
pensate for an error of the LDA and GGA type corre- 
lation functionals that describe effectively only dynam- 
ical correlation coming from electron repulsion at short 
range [15]. Yet the spin-unrestricted KS molecular wave- 
function artifically breaks the symmetry of the disso- 
ciating molecule, displaying unphysical non-zero spin- 
polarization [16] . In fact spurious symmetry breaking has 
remained a much discussed drawback of spin-unrestricted 
DFT (and Hartree-Fock) calculations of molecular prop- 
erties (see e.g. Refs. [17-20]). 

Progress in each of the above noted respects may be 
achieved through orbital-dependent XC functionals ex- 
pressed in terms of the (occupied and unoccupied) KS 
eigenstates, such as the Random Phase Approximation 
(RPA). The RPA functional is the simplest realization of 
the adiabatic-connection fluctuation-dissipation formal- 
ism [21] that defines a broad class of fully nonlocal XC 
functionals. In terms of the "Jacob's ladder" of den- 
sity functional approximations [22], RPA is on the top 
rung, putting it among the most general (but also the 
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most demanding) of present-day approximations. It in- 
cludes the exact KS exchange energy [23] as well as long- 
ranged Coulomb correlations giving rise to dispersion 
forces [6, 24]. This makes the RPA a natural starting 
point to address the above limitations (i) — (Hi) of exist- 
ing functionals in a seamless and consistent way. In this 
paper we show that RPA functional is able to describe 
the strong static correlation in the dissociation of the H2 
bond in a spin-restricted KS formalism, without artificial 
symmetry breaking. The transition from mostly dynami- 
cal to strong static correlation is discussed in terms of the 
adiabatic connection. Analyzing the dissociation energy 
curve of H2, we find that the RPA is accurate around 
the equilibrium bond length and yields, asymptotically, 
the correct dissociation into two H atoms. At interme- 
diate distances the dissociation energy still displays an 
erroneous repulsion. Of course the RPA is just the first 
step in an ongoing systematic quest, with encouraging re- 
sults for small molecules [25, 26] as well as van der Waals 
bonded structures [24, 27, 28]. To achieve satisfactory 
accuracy globally it requires extensions, some of which 
are being examined already [8, 23, 26, 29-31]. 

Last we would like to mention that the XC energy can 
be approximated also as functional of the one-electron 
density matrix, by making an appropriate ansatz for the 
exchange and correlation hole functions. For a particu- 
lar such functional involving both the occupied and the 
unoccupied KS states, Griming et al. [32] recently repro- 
duced the entire dissociation energy curve of H 2 , includ- 
ing proper dissociation. The validity of this approach is 
however still under debate [33]. 

Our paper is organized as follows. In Sec. II we first 
recall the classic problem of stretched H 2 and its implica- 
tions in the context of different density functionals, and 
briefly discuss ways of dealing with it. We also summa- 
rize the RPA equations. In Sec. Ill we examine the XC 
energy of H2 in the framework of the adiabatic connec- 
tion, then, in Sec. IV we present and discuss our results 
for the RPA dissociation energy curve for H2, and com- 
ment on current limitations and ways to overcome them. 
The Sec. V is devoted to results obtained beyond the 
RPA. Section VI summarizes our conclusions. 



II. THE PROBLEM OF STRETCHED H 2 IN 
DFT 

A long-standing problem confronting all single- 
determinant calculations is that of stretched H2, repre- 
sentative for the dissociation of electron pair bonds in 
general [3, 16, 34]. For any bond length R the true 
(interacting) ground state is a singlet [35], with equal 
spin-up and spin-down densities, and the true (noninter- 
acting) KS ground-state corresponds to a single Slater 
determinant \& KS = \® g | made up from the bonding 
<jg molecular orbital. The spin-restricted LDA, GGA, 
and hybrid functionals (as well as Hartree-Fock) correctly 
yield such a ground state, with reasonable total energy, 



around the equilibrium bond length Rq. As a matter 
of fact, the interacting ground state is also mainly of 
|(7 g (7 g | nature around R « Rq. However, as the bond 
length R is increased toward R — > 00, ^ KS no longer re- 
sembles the interacting ground state wavefunction of the 
molecule. Asymptotically the latter assumes the famil- 
iar Heitler-London form and puts precisely one electron 
on each of the two hydrogen atoms, completely suppress- 
ing number fluctuations and describing two free hydro- 
gen atoms, i.e. H •H [16]. By contrast ^ is half 
contaminated by ionic contributions of the form \s a s a \ 
and \sbSb\ (s a and Sb stand for the Is orbitals centered 
on hydrogen A and B respectively), in effect describing 
the stretched H 2 as ± (IT • • • IT ) + \ (H+ • • • H" ) [36] . Op- 
timizing e.g. within spin-restricted GGA, the a g orbital 
does not become a linear combination of the Is hydro- 
genic orbitals but gets much too diffuse, in order to avoid 
the H~ contribution. Hence the dissociation energy of 
stretched H 2 is severely overestimated, its total energy 
lying much above the one of two free hydrogen atoms (as 
illustrated later by our Fig. 5) . On the other hand, a sec- 
ond solution with lower (and reasonable) energy may be 
obtained for bond lengths beyond the so-called Coulson- 
Fishcr point [37, 38] by breaking the spin-symmetry and 
localizing on one hydrogen atom the "spin-up" electron 
and on the other the "spin-down" electron. This is what 
is obtained, in practice, by performing spin-unrestricted 
L(S)DA, GGA, or hybrid functional (as well as unre- 
stricted Hartree-Fock) calculations. 

This situation for approximate XC functionals embod- 
ies the well-known spin symmetry dilemma for dissociat- 
ing H 2 : restricted KS schemes yield the proper symmetry 
adapted ground-state but poor total energies, while unre- 
stricted KS schemes yield very reasonable total energies 
but qualitatively wrong spin densities [16]. 



A. Local, semilocal, and hybrid functionals 

A discussion of this dilemma might appear academic, 
since in most cases one simply calculates the chemically 
bonded system, and subtracts from it the energies of the 
isolated atoms, without ever watching the molecule dis- 
sociate. However, at the origin of the problem is the 
inability of present XC functionals to properly capture 
long ranged left-right correlation that eventually appears 
when a molecule dissociates [39]. Such static (or non- 
dynamical) correlation is present in many real molecules 
even in or near the ground-state geometry. As a mat- 
ter of fact LDA and GGA calculations only mimic it 
through their exchange component while their correla- 
tion component accounts for (short-ranged) dynamical 
correlation only (see e.g. Ref. [15]). In fact their perfor- 
mance depends crucially on the well-known cancellation 
of errors between exchange and correlation [11]. This 
compensation is of course neither perfect nor universal, 
as indicated by the on average overestimated atomization 
energies, particularly for multiply-bonded species like N 2 , 
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but also by often underestimated reaction energy barriers 
(see e.g. Rcf. [40]). 

Hybrid functionals [4, 5] can provide a useful rem- 
edy and are usually more accurate than LDA or GGAs 
alone. Yet they still rely on a similar philosophy: the 
admixing of a certain fixed fraction of the (possibly spin- 
unrestricted) exact exchange energy can be understood 
to improve the description of static correlation while dy- 
namical correlation is still described by (semi-) local LDA 
or GGA functionals. The errors in the exchange and 
correlation energies do not sufficiently cancel in cases 
where either exchange (paradigm: H^) or static correla- 
tion (paradigm: dissociated H 2 ) prevail, so that eventu- 
ally hybrid functionals can become inadequate too [14]. 
For dissociating electron pair bonds, sufficiently nega- 
tive exchange energies are obtained only by resorting 
to spin-unrestricted exchange, i.e. by an unphysical 
breaking of the symmetry of the KS molecular wave- 
function. At which bond length the spin-unrestricted 
solution becomes energetically preferable depends on the 
functional. Hybrid functionals are more prone to break 
symmetry with an increased admixing of exact exchange 
(see Ref. [41]). Artificially symmetry broken solutions 
from spin-unrestricted methods may also yield unphysi- 
cal molecular properties, despite providing a lower energy 
solution and appearing better variationally [18]. We re- 
turn to the concept of hybrid functionals in Sec. Ill A. 



Functionals of occupied and unoccupied KS 
states 



Although the physical origin of the above difficul- 
ties is clear it is far from simple to correct for them, 
while retaining low computational cost. A variety of 
approaches have been applied, mostly along the same 
lines as a traditional restricted Hartree-Fock calculation 
would be corrected, such as spin-unrestricted [38], multi- 
reference [42, 43], or ensemble-referenced Kohn-Sham 
schemes [19, 44]. Although often useful, these suffer from 
similar difficulties as in Hartree-Fock, namely that differ- 
ent approaches work for different situations. 

A more satisfying approach is to develop a more de- 
manding but non-empirical scheme. Functionals involv- 
ing occupied and unoccupied (virtual) KS orbitals of- 
fer this possibility as recently shown by Baerends and 
coworkers [32, 34]. A well-defined starting point is the ex- 
act exchange formalism (EXX) in Kohn-Sham DFT [45] 
which then must be complemented with a compatible cor- 
relation functional, i.e., one that properly respects the 
weakly-correlated, exchange-only limit (e.g. producing 
zero correlation energy in one-electron systems such as 
H^) and accurately interpolates to the strongly- coupled 
regime to be discussed in Sec. Ill A. The RPA XC func- 
tional, discussed next in Sec. II C represents a possible 
first step into this direction. 



C. Functionals from linear response: RPA and 
beyond 

The combination of the random-phase approximation 
(RPA) with DFT was discussed as early as the late 
70's [21], then for the uniform electron gas. The idea 
is to start from the non-interacting KS response function 
and dress it with a (scaled) Coulomb interaction Au co , 
and eventually also with an exchange-correlation kernel 
of time-dependent DFT (TDDFT) [46]. This yields the 
response function of an interacting system, which by way 
of the fluctuation-dissipation theorem gives the corre- 
sponding pair-correlation function and thus the electron- 
electron interaction energy. The XC energy is last ob- 
tained through an integration of the electron-electron in- 
teraction energy along an adiabatic path connecting the 
non-interacting KS system (A = 0) to the interacting one 
(A = 1, see Sec. Ill A). 

Working in the imaginary frequency (iu) domain, the 
basic equations (in Hartree atomic units) are as fol- 
lows. Starting from a KS ground-state with eigenstates 
{4>i a {n}, e i(J [n]}, functionals of the density n, one con- 
structs the KS response 

X u (r,r;iu) = > 

x^(r)^(r)^(r')fe(r') , (1) 

with j ia = 1 for occupied and ji a = for unoccupied 
KS states. The interacting response function \ X at cou- 
pling strength A follows from the Dyson-type screening 
equation 

X \iu) = X °(iu)[l - (Xv cc + /*,(*«)) xV)]^ , (2) 

where v cc — l/|r — r'| is the Coulomb repulsion and 
f£ c (iu) stands for the XC kernel of TDDFT (matrix no- 
tation A =: A(r, r') and AB =: / d 3 r"A(r, r")B(r", r') is 
implied here). The fluctuation-dissipation theorem and 
the coupling strength integration (see Sec. Ill A) finally 
yield an exact expression for the XC energy: 

E xc [n] =-1^a/^_1_ 



•du 



X A (r, r';iu)\ +n(r)5{r-r') 



(3) 



called the adiabatic-connection fluctuation-dissipation 
theorem (ACFDT) for the XC energy. Through approx- 
imations for /xc, fully nonlocal ACFDT XC function- 
als can be generated in practice. In general these will 
be orbital-dependent functionals involving, through x° 
(and eventually /x C ), both occupied and unoccupied KS 
states. The RPA is obtained by setting f£ c = 0, i.e. 
neglecting all XC contributions to screening in Eq. 2. 
Approximating f£ c by the (exact) exchange kernel of 
TDDFT [46] defines the RPA+X functional. In the 
ACFDT formula (Eq. 3), because of the integral over 
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the coupling strength, an approximation to \ °f a given 
order in A, produces an approximation to E xc [n] to the 
next highest order. Thus, to zeroth order, x° inserted in 
Eq. 3 yields the exact exchange energy, E x [n], the first- 
order energy. The RPA yields an approximation to E xc 
that contains all orders of A, but misses contributions 
from second-order onward. In RPA+X on the other hand 
/x C (and thus x A ) are treated exactly up to first-order. 
In turn the RPA+X yields an approximation to E xc [n] 
that is exact up to second-order but, in an approximate 
way, again includes also all high orders of A. Being exact 
to second-order, the RPA+X produces the exact initial 
slope of the adiabatic connection (to be discussed in de- 
tail in Sec. Ill A) as given by second-order Gorling-Levy 
perturbation theory (GL2) [47]. 

While the RPA is expected to be quite accurate for 
(iso-electronic) total energy differences, it tends to un- 
derestimate absolute correlation energies [23]. Indeed, 
the short-range behavior of 

X RPA is 

incorrect, because 

the electrons only respond to the averaged density fluc- 
tuations. In particular, the (spurious) "self-response" of 
an electron to its own contribution to these density fluc- 
tuations is the most important source of error in few elec- 
tron systems. It is responsible for a (spurious) non-zero 
self-correlation energy in one-electron systems like the 
H atom. This deficiency of the RPA may well be cor- 
rected by a local-density functional E^, r-LDA [n] designed 
for the purpose [23] , defining the so-called "RPA+" func- 
tional E*g A+ [n] = El^ A [n} + E s J- hI}A [n} . As expected 
in Refs. [23, 29] and confirmed in Refs. [25, 26], local 
or semilocal short-ranged corrections have rather minor 
effects on molecular dissociation energies, yet nicely cor- 
rect for spurious RPA self-response problem in the H and 
He atoms [26]. For these reasons we focus in this paper 
on the difference of the total energies and XC energies 
between the molecule and the isolated atoms. 

Compared to likewise orbital-dependent Meta- 
GGAs [48-50] and hybrid functionals, the RPA involves 
both occupied and unoccupied KS states. Its computa- 
tional cost is quite severe (a factor 10 2 . . . 10 3 compared 
to a GGA approach), but its promise is to tackle the 
nonlocality of exchange and correlation on an equal 
footing, leading to a seamless description from the 
chemically-bonded to the dissociated regime, including 
van der Waals interactions not accounted for in GGAs 
or hybrids. 

Recent calculations for small molecules by Furche [25] 
found that the RPA describes the chemical bonds with 
similar but not better accuracy as modern GGAs. Al- 
though this appears to be disappointing, we note that 
the RPA is in fact accurate for H2 and the difficult case 
of Be2 [26] where LDA and GGA fail. Also there has been 
significant progress in building XC functionals on top of 
the RPA that include the dispersion forces between lay- 
ered solid-state structures such as jellium slab models [24] 
and graphite [27], and between atoms or molecules [28]. 
We further point out that the RPA is just a ready real- 
ization of a much broader class of functionals based on 



the adiabatic-conncction fluctuation-dissipation theorem 
which offers various options for improvements [26, 29- 
31]. Moreover the RPA is amenable to extensions de- 
rived in many-body Green function theory which itself 
is being very actively explored for total energy calcula- 
tions [10, 51, 52], including H 2 [53-55]. 

In the present study we do not perform self-consistent 
RPA calculations, these are computationally too costly 
at present (although formally perfectly feasible, see 
Refs. [56, 57]). Instead we evaluate Eqs. 1 - 3 a posteriori 
with the KS eigenstates taken from EXX calculations as 
explained further in Sec. HID. Thanks to the variational 
principle for DFT total energies we expect the resulting 
RPA total energies to be tight upper bounds to the self- 
consistent RPA results. Indeed we found that our results 
for the RPA total energies remained virtually unchanged 
when we used the LDA [58] or PBE GGA [59] instead of 
the EXX for calculating the initial KS states. 

III. ANALYSIS OF THE H 2 DISSOCIATION 

In this section and in the next one, we show that the 
RPA offers a promising handle on the problem of dis- 
sociating electron pair bonds as exemplified by H 2 . In 
particular, we give the correct prescription for applying 
the scheme during dissociation and show that its remain- 
ing errors may well be correctable. Using it, we also 
provide the first accurate calculations of adiabatic con- 
nection curves as a system passes through its Coulson- 
Fishcr point where the KS ground-state for standard XC 
functionals bifurcates in a spin-symmetry adapted and a 
lower energy but symmetry-broken solution. 

A. Adiabatic connection 

To understand in detail how DFT handles static cor- 
relation, one invokes the adiabatic connection [21, 38], 
already briefly mentioned in Sec. II C. One imagines al- 
tering the strength of the electron-electron repulsion by 
multiplying it by a constant A, which varies between 
and 1. At the same time, the one-body potential is al- 
tered, i.e., made a function of A, so as to keep the electron 
density fixed. This is a way of continuously connecting 
the noninteracting KS system (A ~ 0) to the interacting 
physical system (A = 1). More importantly, by virtue of 
the Hellmann-Feynman theorem, one can write the XC 
energy as an integral over purely potential energy: 

E xc [n] = f d\ U xc [n](X) (4) 
Jo 

where U xc [n}(\) = (^ x [n]\v ee \^ x [n]) - U[n}. Here * A [n] 
is the ground-state wavefunction at coupling strength A, 
v cc is the Coulomb repulsion, and U[n] is the Hartree 
energy. The integrand U xc [n](X) represents the potential 
energy contribution to XC and makes up the adiabatic 
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connection curve. At the A = end it corresponds to the 
exact Kohn-Sham exchange energy E x [23], 



U xc [n}(0) = E x [n] 



(5) 



and has a (negative) initial slope, given by the correlation 
energy of second-order Gorling-Levy perturbation theory 
(GL2) [47] 



K a [n](X) = —U xo [n](X) 



= 2E^[n] . (6) 



A=0 



All A-dependence rests in the correlation contribution 

U c [n}(X) = U xc [n}(X)-E x [n] . (7) 

At A = 1, U xc [n](X) describes the XC potential energy of 
the physical system, U xc [n] =: U xc [n](l), and, similarly, 
the correlation potential energy, f/ c [n] =: £/ c [n](l). Since 
correlation reduces the electron-electron repulsion, the 
different energies are always ordered as 



E x [n] > E xc [n] > U xc [n] 



(8) 



The A-dependence for ACFDT type functionals (see 
Eq. 3), appearing through the response function % A , is 
given by 



U c [n](X) = - 



Tr[ Wco { X A ( lU )-x°M}] 



(9) 

where Tr[A] =: J d 3 rA(r,r) . The full A-curve may then 
be calculated thanks to Eq. 7. For the LDA or GGA type 
functionals the A-dependence can be readily calculated 
from the exact relation [60] 



d 



U* c [n}(X) = — (A 2 £ xc [n 1/A ]) 



(10) 



where the XC energy functional is evaluated at the scaled 
density ni/\{v) =: n(r/A)/A 3 . Hence analysis by adi- 
abatic decomposition allows investigation into how well 
the different functionals perform, and why [11, 61]. 

Below we consider molecular dissociation energies, i.e. 
the difference between the molecular and atomic total 
energies, 

AE = E tot [molecule] — E tot [atoms] , (11) 

and analyze the exchange-correlation contributions by 
means of the analogously defined differences AE XCl AE X , 
and AU XC (X). 

A useful measure of the correlation strength is given 
in Rcf. [62]. Define the parameter b by 



E xc [n] = bE x [n] + {l-b)U xc [n] 



(12) 



A simple interpretation of b is given by the following 
geometrical construction: if the adiabatic curve were a 
horizontal line of value -E x [n] running from to 6, and 
then dropped discontinuously to another horizontal line 
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FIG. 1: Adiabatic connection for H2 at bond length R = 
1.4 bohr within the RPA (solid line) and the GGA (dot- 
dashed line). The GL2 curve (dashed) corresponds to the 
slope of the exact curve at A = 0. Shown is the difference be- 
tween H2 and two free H atoms, evaluated on self-consistent 
EXX densities. The "exact" curve is an interpolation [65] 
based on accurate values of AE X , AE XC , and AU XC from a 
configuration interaction calculation [63, 64]. 



of value t/jccM running from b to 1, then b is that value 
of A that yields the correct E xc [n]. Thus, in the high 
density limit, where the adiabatic connection curve is a 
straight line, b is exactly 1/2. On the other hand, for 
strong static correlation, in which the adiabatic connec- 
tion curve drops rapidly to its final value, b is close to 
zero. One can also show [62] 



b = 



T c [n] 

\UrM\ 



(13) 



where T c is the kinetic portion of the correlation energy 

T c [n) =E c [n]-U c [n} . (14) 

Thus small b indicates that the correlation is indeed 
static, i.e., has a smaller fraction of kinetic to potential 
energy. 

For the atoms and most chemically bonded systems, 
the adiabatic connection curve is rather non-descript, ly- 
ing close to a straight line. This is illustrated by H2 
at the equilibrium bond length in Fig. 1, where we plot 
AU X ^ A (X). The area enclosed by the adiabatic connec- 
tion curves represents the XC contribution to the dis- 
sociation energy, AE XC . As can be seen from Tab. I 
the RPA dissociation energy of H 2 is excellent agreement 
with the exact value. Noting also the good agreement of 
the endpoints of our RPA curve in Fig. 1 with accurate 
data from configuration interaction calculations [63, 64], 
listed in Tab. I, this implies that the RPA curve lies very 
close to the true adiabatic connection curve. The straight 
line corresponds to GL2 theory, indicating the 6=1/2 
high density limit and the initial slope of the true curve 
given by Eq. 6. For A > 0, Af/° L2 (A) lies below the 
true curve, overestimating the absolute correlation en- 
ergy. Indeed, calculating ^ RPA +x to ^ rst or d cr in A we 
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TABLE I: Adiabatic decomposition of the dissociation energy 
AE of H 2 at bond length R — 1.4 bohr, evaluated on self- 
consistent EXX densities. Shown are the differences between 
the molecule and two free H atoms for the coupling strength 
integrand AUxc(^) and related quantities, as explained in the 
text. All values are in eV. 





AE 


AE XC 


AE X 


AU C 


AT C 
|Al/ c 


At/ xc (0) 


PBE GGA 


-4.54 


-1.91 


-1.04 


-1.52 


0.431 


-2.42 


RPA 


-4.73 


-2.10 


-0.99 


-1.95 


0.427 


-2.97 


exact 


-4.74" 


-2.04 6 


-0.98 6 


-1.93 6 


0.450 


-2.84 



"Reference [35]. 

'From Rcfs. [63, 64]. In these works, E x , E x c, and Tq were calcu- 
lated on the H2 density obtained from a configuration interaction 
calculation which yielded AE = —4.68 eV. Using these data we 
evaluated AUc from Eq. 14. 



obtain AE GL2 = -5.04 eV, about 0.3 eV below the true 
dissociation energy. This indicates that the second-order 
perturbative treatment is qualitatively but not quantita- 
tively accurate for H 2 . Regarding the PBE GGA, Fig. 1 
and Tab. I show that it is accurate for the exchange en- 
ergy (A = 0) but underestimates the absolute correlation 
potential energy (A = 1). In turn the PBE GGA also un- 
derstimates the absolute XC energy and the dissociation 
energy of H2 . 

Calculations of the exact adiabatic connection, us- 
ing the accurate ground-state wavefunctions for all A, 
have so far not been attempted to our knowledge for 
molecules, while a few such curves based on accurate 
ground-state densities have been reported for atoms [66- 
68], bulk Si [69], and model systems [70]. We remark 
that £U xc [n}(\)\ A=cb i- c - the GL2 correlation energy, is 
also a key ingredient in a recent coupling strength inter- 
polation of the adiabatic connection by Seidl et al. [71] 
which performs with similar accuracy as modern hybrid 
functionals for molecular dissociation energies. For H2 
close to the equilibrium bond length, the corresponding 
dissociation energy and adiabatic connection curve are in 
good agreement with our RPA results [72]. 



B. H2 symmetry dilemma 

We now discuss the stretching of H2 as a paradigm 
of the difficulties that single-determinant methods have 
with dissociation. The nearly straight line behavior dra- 
matically changes when the bond is stretched to R — > 00. 
Asymptotically, the proper molecular wavefunction for 
any A > (i.e. regardless of the interaction strength) is 
the bonding linear combination of the Is orbitals of the 
two H atoms. For the H 2 "super- molecule" the exchange 
energy therefore has the same value as in a hydrogen 
atom, i.e. £ X [H---H] E X [E] = -U[n K ]- Conse- 
quently also U C [H - ■ ■ H'](A) = — £7[nn] f° r an Y A > 0, 
since the dissociated H2 molecule must have the same 
total energy as two H atoms. Figure 2 shows the corre- 
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FIG. 2: Exact adiabatic connection for dissociated H2 at bond 
length R — ► 00, shown as the difference with respect to two 
free H atoms. For the exact KS determinant the curve starts 
at A = with the negative of the exact exchange energy of 
a single H atom, — E X (H) ~ 8.5 eV as explained in the text. 
The negative of the shaded area represents the correlation 
energy of the H2 molecule and equals E X (H). 



sponding exact adiabatic connection AU XC [R'- ■ -H](A). 
The immediate drop of the adiabatic connection curve at 
A = is characteristic for a system with strong static cor- 
relation: here 6 = exactly. The position of one electron 
entirely determines the position of the other electron: the 
two electrons in infinitely separated H2 must sit on the 
two different nuclei but never on the same, as spuriously 
allowed by the single KS determinant. This type of cor- 
relation is often refered to as left-right or static corre- 
lation, where the true many-electron wavefunction takes 
on multi-determinant character [3]. Put differently [34], 
in the concept of the XC hole, the exchange hole of H2 is 
spatially completely delocalized over both nuclei. How- 
ever the XC true hole is always centered about the refer- 
ence electron. This means that the correlation hole must 
be long-ranged to yield the proper hydrogen-like hole on 
one nucleus and the needed zero total XC hole on the 
opposite nucleus. By contrast LDA or GGA correlation 
holes, derived essentially from the uniform electron gas, 
are always short-ranged, and hence cannot cancel the ex- 
act exchange hole far away. Breaking inversion symme- 
try, L(S)DA and spin-dependent GGA on the other hand 
(like unrestricted Hartree-Fock) already yield the spin-up 
and -down exchange holes of separate hydrogen atoms, 
opposite spin electrons sitting on different nuclei, and 
thus mimic the static correlation: unrestricted Hartree- 
Fock or exact KS exchange indeed yield two hydrogen 
atoms as the dissociation products, and produce curves 
similar to Fig. 2. 



C. Behavior of hybrid functionals 

Does any of this matter for real systems, especially if 
we never ask them to dissociate? The answer is emphati- 
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FIG. 3: Sketch of AU XC (\) for N 2 -> 2N. (The EXX, 
PBE, and PBEO hybrid curves are calculated from the self- 
consistent PBE densities and orbitals. The "exact" curve is an 
interpolation between A_E X and A(7 X c E , following Ref. [61].) 

cally yes. Currently our most accurate popular function- 
al are hybrids, which mix a fraction ao s=a 1/4 of exact 
exchange with GGA. To understand why this improves 
the accuracy over GGAs [11], one first notes that usually 
the latter work best near A = 1 (H2 being an exception 
to this rule) [11]. But they worsen as A — > 0, especially 
when there is static correlation in the system. The reason 
is twofold and illustrated in Fig. 3 for N2 (where static 
correlation is expected to be most pronounced in the ir 
bonds). On the one hand GGA exchange energies tend 
to be too negative for molecules, leading to a too low 
starting point for the adiabatic connection at A = 0. On 
the other hand the GGA correlation potential energy is 
accurate toward A = 1, but too flat a function other- 
wise. As a consequence Al/ XC (A) in the GGA lies too 
low, leading to an overestimate of the magnitude of the 
dissociation energy. Strong static correlation shows up 
as a rapid drop of J7 XC (A) from A = toward A = 1. 
Note that the free N atoms contain little static corre- 
lation, that is, their U XC (X) lie close to a straight line 
and can be accurately obtained from GL2 perturbation 
theory. Recognizing that GGAs (and also the LDA) are 
typically least accurate at the A = end, hybrid func- 
tionals can improve over the GGA adiabatic connection 
by mixing in the exact exchange energy, i.e. part of the 
exact limit U xc (0) = E x . This has led to hybrid XC 
functionals of the generic form 

E^ c h ^E^ A +a (E x -E^ GA ) • (15) 

Becke has determined an emprirical value of ao s=s 0.28 
[4]. The non-empirical estimate of ao = 1/4 has been 
given, based on the observation that ?7 X c(A) is formally 
provided by DFT perturbation theory and fourth order is 
adequate for typical molecules, and defines the so-called 
PBEO hybrid [5], shown for N2 in Fig. 3. Hybrid function- 
als like the PBEO on the average yield molecular dissoci- 
ation energies close to chemical accuracy [73] , improving 
over GGA. Pictorially, the admixing of exact exchange 



shifts the GGA curve a fraction ao toward the true start- 
ing point (A = 0) of the correct curve. Since the true 
adiabatic connection curve drops sharply from A = 0, a 
fraction ao < 1/2 produces a hybrid curve whose area is 
about that of the true curve, as in Fig. 3. 

From the adiabatic connection we can also understand 
when (semi-)local and hybrid functionals become inap- 
propriate. If the XC energy is dominated by exchange 
and there is little correlation energy, the adiabatic con- 
nection curve drops little from its A — value. For one- 
electron systems, such as in , there is even zero cor- 
relation and the adiabatic connection curve stays com- 
pletely flat. Then one ought to use ao w 1, i.e. 100% 
exact exchange. Indeed GGA as well as hybrid func- 
tionals perform poorly for one-electron bonds such as in 
radical ions [12, 13]. On the other hand if static corre- 
lation is important, the adiabatic connection will drop 
rapidly from its value at A = 0, given by the exact ex- 
change energy. This drop is missed by GGA correlation 
functionals, which then underestimate the correlation en- 
ergy. This error may be compensated by using approxi- 
mate GGA exchange functionals since these tend to over- 
estimate the absolute exchange energy in molecules and 
the exchange contribution to dissociation energies [11]. 
Hybrids with ao > improve the exchange energies over 
GGA alone but eventually also reduce the error cancella- 
tion due to GGA exchange. An early indication for this 
was the observation that combining 100% exact exchange 
with GGA correlation completely fails for molecular dis- 
sociation energies [74]. In fact no exact exchange mixing 
(ao ~ 0) can sometimes be the better choice as demon- 
strated in Ref. [14]. We stress that to describe strong 
static correlation effects in bond dissociation, any of these 
functionals must be used in spin-unrestricted form in or- 
der to attain exchange energies that are negative enough 
to compensate for the lack in the respective approxima- 
tion for correlation. Clearly, to describe the transition 
from exchange dominated to static correlation dominated 
XC regimes, a more flexible, i.e. "system sensitive" in- 
terpolation of the adiabatic connection than provided by 
present hybrid functionals is needed [62]. 



D. Beyond the Coulson-Fisher point 

A key concept in this paper is that DFT within the 
RPA allows correct dissociation of molecules. However, 
given our present inability to perform self-consistent RPA 
calculations, the demonstration of this fact becomes quite 
subtle. While the (restricted) EXX solution is an ade- 
quate starting point for the RPA around the equilibrium 
bond length, it is totally inadequate beyond the Coulson- 
Fischer point (R > 2.5 bohr for H 2 treated in EXX), 
where ambiguity arises in a single-determinant calcula- 
tion. In the words of the symmetry dilemma, should one 
use the unrestricted solution, which has a pretty good en- 
ergy but totally incorrect spin density, or the restricted 
solution, which has the correct symmetry but poor ener- 
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FIG. 4: Same as Fig. 1, but for R = 5 bohr, i.e. beyond the 
Coulson-Fisher point. The RPA results are based on the total 
density of a unrestricted EXX KS calculation. Also shown 
are the adiabatic connections for the PBE GGA applied in 
the restricted KS formalism (RKS), yielding poor energetics, 
and in the unrestricted KS formalism (UKS) , yielding better 
energetics but artifically breaking inversion symmetry. 

getics (as seen in Fig. 5) ? 

The answer has been given in Ref. [16]. One must 
use the best estimate for the correct ground-state den- 
sity that is available. As argued there, the unrestricted 
solution yields the best approximate DFT density, but its 
spin density is not to be believed. We therefore take the 
total density from our unrestricted EXX KS calculation, 
and treat it as a spin-singlet. This becomes our input 
density to our RPA calculation. Recall that this density 
becomes exact in the limit of R — ► oo, where it corre- 
sponds to two separate hydrogenic densities, in contrast 
to a restricted scheme. Inverting the KS equation [75] for 
n(r) = n^ xx (r) + nf xx (r), we obtain the KS potential 
yielding that density and the respective KS eigenstates. 

IV. H 2 DISSOCIATION WITHIN THE RPA 

As can be seen from Fig. 4, the RPA adiabatic con- 
nection of H 2 beyond the Coulson-Fisher point becomes 
strongly bent downward, capturing the onsetting crucial 
contribution of static correlation related to the multi- 
determinant nature of the interacting many-electron 
wave function. The latter is not captured in the PBE 
GGA correlation functional which significantly underes- 
timates the magnitude of the correlation energy, as seen 
from Tab. II and evident from the PBE GGA (RKS) 
curve in Fig. 4. Switching to the unrestricted PBE (and 
PBEO) schemes, their exchange component simulates the 
missing static correlation while the correlation compo- 
nent is much too small in magnitude as Tab. II shows. 

Correspondingly, unrestricted PBE (like PBEO) even- 
tually give dissociation energies that are in quite good 
agreement with the exact value, as listed in Tab. III. Wc 
stress that this is a result of error cancellation between 
(unrestricted) exchange and correlation: the PBE adi- 



TABLE II: Adiabatic decomposition of the dissociation en- 
ergy AE of H2 at bond length R = 5 bohr, as shown in 
Fig. 4. The RPA functional is evaluated on the total den- 
sity from a spin-unrestricted EXX calculation as described in 
Sec. HID. Also shown are results for the PBE GGA, evalu- 
ated as a spin-restricted (RKS) and a spin-unrestricted (UKS) 
functional using the EXX spin-densities. All values are in eV. 





AE 


AExc 




AU C 


AT C 


A(7ic(0) 


PBE (RKS) 


2.20 


2.30 


2.81 


-0.90 


0.43 


-1.23 


PBE (UKS) 


-0.06 


0.040 


0.044 


-0.008 


0.50 


0.084 


RPA 


0.54 


1.33 


5.72 


-5.06 


0.13 


-111 


exact 


-0.10" 


0.82 6 


5.85 6 


-5.61 6 


0.10 


-56.8 



"Reference [35]. 

'From Refs. [63, 64], see also Tab. I. 



abatic connection curves is qualitatively clearly wrong, 
especially at A = 0, where its A.E X is much too small and 
its slope turns out even slightly positive (though this is 
not visible on the scale of Fig. 4). No such error cancel- 
lation occurs within the RPA. The onset of strong static 
correlation is reflected in the low value of our correlation 
strength parameter b reported in Tab. II for the RPA and 
exact curves, but not for (semi-) local density functional 
approximations. 

Although qualitatively correct the RPA adiabatic con- 
nection curve is still deficient, as can be appreciated by 
comparing to the exact curve in Fig. 4 and data in Tab. II. 
AIT? PA (A) does not drop deep enough with A, despite its 
too steep initial slope. Correspondingly the RPA yields 
a too positive molecular correlation energy and produces 
an artificial barrier for dissociation as seen in Tab. III. 
This is further evidenced in the full RPA dissociation 
curve of Fig. 5. Indeed, while the RPA performs accu- 
rately around equilibrium R and again at larger R, it 
shows an unphysical bump at intermediate bond length 
R. The origin of this bump will be further discussed be- 
low. Nonetheless the asymptotic behavior (R — » 00) of 
the RPA is correct, as can be understood from a model 



TABLE III: Dissociation energy of H2 at bond length 7?, cal- 
culated with different XC functionals as indicated in Fig. 5. 
Given are the energies from unrestricted KS calculations, ex- 
cept for RPA and RPA+X as explained in the text. All values 
are in eV. 
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EXX 


-3.62 


-0.47 


-0.02 


0.00 


PBE GGA 


-4.53 


-1.27 


-0.03 


0.00 


PBEO hybrid 


-4.52 


-1.07 


-0.03 


0.00 


RPA 


-4.73 


-1.44 


+0.54 


+0.20 


RPA+X 


-4.86 


-1.45 


+0.34 


-0.25 


exact a 


-4.75 


-1.56 
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0.00 



"Reference [35]. 



9 




bond length (bohr) 

FIG. 5: Dissociation energy curve of H2. The upper panel 
compares results for EXX, the PBE GGA, and the PBEO 
hybrid functionals calculated in the self-consistent restricted 
(RKS) and unrestricted (UKS) KS formalism, with exact 
data from Ref. [35]. The lower panel compares the RPA and 
RPA+X curves, calculated with total densities n UKS obtained 
from unrestricted EXX KS calculations, with exact data. Also 
shown is the RPA curve calculated with densities n RKS from 
restricted EXX KS calculations. 



RPA calculation using only the HOMO and LUMO EXX- 
KS states. As shown in the Appendix, this model RPA 
calculation yields the exact correlation (and total) en- 
ergy for R — * 00, and produces precisely the exact adi- 
abatic connection of Fig. 2. Including the higher lying 
KS states, the RPA builds up spurious self-correlation in 
both the H atoms and the H2 super-molecule, which how- 
ever cancels out in the dissociation energy. This cancel- 
lation is indeed also reflected in the (identical) estimates 
of the local-density corrections (RPA+) in the atom and 
in the infinitely stretched molecule. 

Figure 5 also demonstrates that for proper dissocia- 
tion it is crucial to work with qualitatively correct densi- 
ties, i.e., to start from a KS potential that takes into ac- 
count the essential left-right correlation. Describing the 
system by the restricted KS formalism using the EXX, 
PBE GGA or PBEO hybrid functionals leads to much 
too high total energies for the dissociating bond because 
the self-consistent EXX density is qualitatively wrong as 
explained in Sec. II. While the PBE GGA and PBEO 
include correlation, improving over EXX, it is obvious 



from Fig. 5 that only the unrestricted KS formalism leads 
to reasonably accurate dissociation energies for larger R. 
The fact that the PBEO curve raises too quickly above the 
true curve indicates that the error cancellation between 
exchange and correlation is significantly incomplete, as 
explained in Sec. Ill A. For H2, PBE GGA clearly works 
better, and the lower energy solution appears only be- 
yond R s=s 4 bohr, i.e. for a markedly larger bond length 
than in both the EXX and the PBEO hybrid. In Ta- 
ble III we list the dissociation energies for the different 
functionals at various bond lengths. 

The RPA curve is accurate around the H 2 equilib- 
rium bond length and approaches the dissociation limit 
for large R. The success of the RPA lies in the fact 
that it does so (properly) as a functional of a singlet 
density only rather than spin densities as in the tra- 
ditional approximations for XC. Of course the density 
must ultimately come from a self-consistent KS calcu- 
lation and potential, whereas we have approximated it 
non-selfconsistcntly. Our findings confirm that proper 
densities require accurate approximations also for the XC 
potential (the functional derivative of E xc [n]), as argued 
recently by Baerends [34]. In agreement with Ref. [34], 
our results also show that unoccupied KS states (and in 
particular the LUMO state) must be included in E xc [n] 
in order to attain the correct dissociation limit. 

An obvious shortcoming of the RPA compared to the 
unrestricted PBE GGA and PBEO schemes is that its 
dissociation curve displays an unphysical bump for in- 
termediate R as seen in Fig. 5 and Tab. III. A similar 
behavior has been found for N 2 [25] and Be 2 [26, 76]. We 
believe that these deficiencies stem from the RPA itself, 
rather than our lack of self-consistency. In particular we 
have obtained essentially the same curves when we used 
densities derived from unrestricted LDA and GGA cal- 
culations. A more conclusive answer must either start 
from more accurate densities or await self-consistent cal- 
culations. Using an ansatz for the XC energy functional 
in terms of natural orbitals in a (approximately self- 
consistent) singlet KS calculation, Baerends and cowork- 
ers [32, 34] recently reported a dissociation curve of H 2 
in very good agreement with the true curve for all R. 
The fact that for intermediate R their curve is distinctly 
more accurate than our non-self-consistent RPA result 
clearly calls for further analysis of both approaches, in- 
cluding possible error cancellation between different com- 
ponents of the total energy. This is beyond the scope of 
our present study. 



V. EXTENSIONS BEYOND THE RPA 

So far we have shown that the RPA gives a qualita- 
tively correct account of the (differential) adiabatic con- 
nection of H 2 , in contrast to semilocal or hybrid func- 
tionals, but needs to be improved for moderately large 
bond lengths. We now discuss possible extensions of the 
RPA as an ACFDT XC functional. 
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A known deficiency of the RPA is the spurious self- 
correlation in the absolute correlation energies, for the H 
atom as well as the H2 molecule. For the one-electron 
H atom, self-correlation is eliminated by including the 
exact exchange kernel (/£ = — Au co ) in the screening of 
the electronic Coulomb interaction, Eq. 2. For the two- 
electron H 2 the exact exchange kernel (/* = — ^v cc [46]), 
eliminates self-correlation to second order in A yet not to 
higher order. Using this RPA+X functional we obtain 
the dissociation energy of H 2 (R = 1.4 bohr) as —4.86 eV, 
about 0.1 eV below the true and our RPA values. Having 
estimated our computational accuracy at the same order, 
we feel cautious about the significance of the RPA+X re- 
sult. On the other hand we find that the dissociation 
energy curve beyond the Coulson-Fisher point is not at 
all improved. To the contrary, while the RPA+X curve 
follows the RPA curve up to w 4 bohr, it drops below 
zero and approaches a negative constant for larger R. 
We regard this as a size-consistency problem in that the 
self-correlation error is eliminated in the H atoms, but 
reappears in the far stretched H 2 . This is further cor- 
roborated in Appendix B. Our finding clearly suggests 
that in addition to exchange also correlation contribu- 
tions need to be included in the XC kernel. A ready 
way to do so would be to employ the well-known adia- 
batic LDA kernel [77] or an energy-optimized adiabatic 
XC kernel [30]. 

As mentioned in Sec. II C the incorrect RPA short- 
ranged correlation, including self-correlation, may be cor- 
rected through a specially designed (semi-)local-dcnsity 
functional (RPA+). Although this improves upon the too 
negative RPA correlation energies as shown in Ref. [26], 
it again does not correct the deficiencies we observe in 
the RPA dissociation energy curve of H 2 . 

The limitation of RPA, as in any adiabatic treatment 
of the interacting linear response, might be that it treats 
only single excitations [78] and thus cannot take into 
account contributions to density fluctuations that cor- 
respond to doubly excited determinants and eventually 
also contribute to the correlation energy. As is well 
known, doubly excited determinants have larger weights 
in the asymptotic interacting wavefunctions of H 2 for 
both ground and excited states. For R — > 00 this is evi- 
dent from the corresponding exact Heitler-London wave- 
functions (see e.g. Refs. [79] and [80]; for instance, the 
lowest 1£+ and 1S+ singlets are the symmetric and an- 
tisymmetric linear combinations, respectively, of the de- 
terminants Icr^CTgl and \a u a u \ made up of the HOMO 
and LUMO). Double excitations imply additional poles 
in the (real frequency) interacting response and thus a 
strongly frequency-dependent XC kernel. Any such pole 
contributes to the spectral decomposition of the pair- 
correlation function or XC hole. While there is limited 
progress concerning the calculation of the excitation ener- 
gies for certain doubly excited states within TDDFT [81], 
further analysis (and development) of the spatial and fre- 
quency dependence of such XC kernels is needed for ap- 
plications in ACDFT XC functionals. 



VI. SUMMARY 

The central message of this paper is that the Random 
Phase Approximation (RPA) resolves the long-standing 
symmetry dilemma encountered in approximate density 
functional theory when breaking the H 2 electron pair 
bond. We showed that the RPA produces the correct 
dissociation limit from a proper singlet KS density, 
without the need for artifical symmetry-breaking as in 
unrestricted Kohn-Sham theory for traditional local- or 
generalized-gradient functionals and hybrid exchange- 
correlation (XC) functionals. By analyzing the adiabatic 
connection, we showed that the RPA captures correctly 
the strong static (left-right) correlation that arises when 
the pair bond breaks. Local and gradient-corrected func- 
tionals make serious errors here, and even hybrids, which 
mimic this effect at equilibrium bond lengths, cannot 
account for the extreme static correlation in the disso- 
ciation limit. As the RPA yields an orbital-dependent 
XC functional, our results demonstrate the importance 
of including unoccupied KS states, in particular the 
LUMO state. We also showed that it is crucial to work 
with an accurate density, which we have constructed 
approximately in this study, and that could be improved 
by applying the RPA self-consistently. We further found 
the RPA dissociation curve to agree well with exact data 
near the equilibrium bond length of H 2 . When the bond 
is stretched, it tends to the correct limit, unlike all the 
common restricted Kohn-Sham approaches. Noting that 
the RPA still leads to an unphysical repulsion of the 
hydrogens for intermediate bond lengths, we adressed 
inherent limitations and possible extensions of the RPA. 
Seen as a first step to realize fully nonlocal XC function- 
als by the adiabatic-connection fluctuation-dissipation 
formalism, we believe that the RPA provides a sound 
basis for quantitative refinements. Our study highlights 
H 2 clS cl SI gnificant benchmark system for assessing 
future progress beyond the RPA. 
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Appendix 

In part A of this Appendix we describe our computa- 
tional method for evaluating the RPA XC energy. In 
part B we corroborate that the RPA total energy of 
H2 at asymptotically large bond distances from a spin- 
restricted Kohn-Sham groundstate (i) becomes indeed 
equivalent to the total energy of two free H atoms as cal- 
culated in a spin-unrestricted formalism, in contrast to 
the case of the RPA+X kernel, (ii) for sufficiently large 
separations R includes the expected —Cf PA /R Q van der 
Waals attraction. 



A. Computational Method 

We have implemented the RPA functional in a pseu- 
dopotential plane- wave framework [82], handling the re- 
sponse functions, Eqs. 1 and 2, in their reciprocal space 
representation. Gauss-Legendre quadrature rules are 
used for the A and iu integrations. For our study of 
H 2 the — 1/r attraction is replaced by a highly accurate 
norm-conserving pseudopotential [83]. The latter yields 
practically the exact energy of the H atom and disso- 
ciation energies to within 0.1 mHa when compared to 
full-potential results, for LDA, GGA, or EXX calcula- 
tions. We place the H2 molecule in a fee supercell of 
21 bohr side length. For the initial KS calculation we 
use a plane- wave cutoff energy of 30 Ha, and 12 Ha for 
the response functions. In the KS response we include 
unoccupied states up to 2.5 Ha explicitly and treat the 
higher ones through a closure relation. For the frequency 
integration we employ 12 supports for R < 2 bohr and 
up to 54 for larger R, concomitant with the closing of 
the HOMO-LUMO gap. For the coupling strength inte- 
gration we use 4 to 11 supports to capture the stronger 
curvature of U£ c with increasing R. From convergence 
tests we estimate that our total and dissociation energies 
are converged to well within 0.1 eV. Indeed our value 
for the H 2 dissociation energy in RPA, —4.73 eV, is in 
excellent agreement with previous work [25, 26]. 

B. RPA XC for H2 at large bond lengths 

In this section we show analytically that within the 
RPA lim^ 00 S t ot(H 2 ) = 2£ tot (H), i.e. that the total 
energies of the infinitely stretched, spin-compensated H 2 
molecule and of two separate spin-polarized H atoms are 
identical. We also show that this does not hold for the ex- 
act exchange kernel (RPA+X approximation). We will 
examine the density response of the H 2 molecule using 



the particle-hole formulation (in a product basis of the 
KS states) of time-dependent DFT [78], which is equiva- 
lent to the matrix formalism used in Sec. II C but more 
convenient for formal analysis. We emphasize that our 
analysis holds for the RPA as an adiabatic approximation 
only, and does not address the effects of double excita- 
tions that require a non-adiabatic (frequency-dependent) 
correlation kernel as argued in Sec. V. 

For the spin-compensated molecule the response func- 
tion for some coupling strength A can be written as 

^(ry;«)=E«iWra)Oi'M . as) 

n n\ J 

where ft n (X) are the transition frequencies and are 
the associated amplitudes. The f2„(A) are the positive 
square roots of the eigenvalues of the matrix 

M ijM (A) = uljSikSjt + 4Vw«/™ fc °(A)v^M (17) 

involving all possible KS excitations from an occupied 
KS state 4n to an unoccupied KS state <f>j , with re- 
spective KS transition frequencies = Cj — e^. Here 
m°lW = /d 3 rd 3 r>,(r)^(r)/ H \ c (r,r')^(r')0/(r') de- 
notes the matrix elements of the Hartree and XC kernel, 
and it is assumed that one works in the adiabatic approxi- 
mation, i.e. that the XC kernel is frequency independent. 
From the eigenvectors n of M^^A), one obtains the 
spectral components 

q£. n = yW^W u* n (is) 

of the amplitudes Qn( r ) = 

E occ. E unocc. QA. n # (r) ^. (r) . T n thc RPA, / H \ c = Xv cc . 

Consider now the H 2 molecule at large R. For any 
finite number TV of (bound) KS states, R can be cho- 
sen large enough such that the molecular orbitals can 
be approximated by the bonding and antibonding linear 
combinations of the atomic orbitals a,(r) and 6j(r) of the 
H atoms A and B respectively: 

&±(r)~ (2T2^)- 1/2 { ai (r)±6,(r)} , (19) 

where Si is the overlap integral. Similarly, we approxi- 
mate the respective eigenenergies Ei± ~ a by the atomic 
eigenenergies q , except for the HOMO and LUMO ener- 
gies Eo± whose gap we write as E g = Eq- — Eq + . 

In the (imaginary frequency) KS response x°[H 2 ] of the 
molecule, we first split off the HOMO-LUMO transition 
which is well separated from thc higher ones, and define: 

5 x O[H 2 ](r,r',i«) = -4£7 fl g^ , (20) 

where q(r) = o+ (r)0o-(r) - (2^T^)- 1 (\a (r)\ 2 - 
|6o(r)| 2 ). For R — ► 00, the remainder x°[H 2 ] =: x°[H 2 ] — 
5x°[H 2 ] can readily be shown to split (up to exponentially 
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decreasing corrections) into atomic contributions x° [Ha] 
and x°[H B ], where: 

o [TT , • , v^,, ^(^-(r^'Mr') 

(21) 

with a similar definition for x°[Hb]- X°[Ha] and x°[Hb] 
are formally equivalent to the KS reponse of a free spin- 
polarized H atom (see however the discussion of RPA+X 
below). Hence we asymptotically get: 

X°[H 2 ] = x°[H A ] + x°[H B ] + ^X°[H 2 ] + 0(cxp.) . (22) 

As for the asymptotic RPA response (R — > oo), inspec- 
tion of Eq. 17 shows that the lowest eigenvalue Oq(A) (i.e. 
the singlet excitation energy) exponentially tends to zero 
and is well separated from the others. Indeed, 

fig(A) = E 2 g + AXK E g + X 2 P(X)E g + O {{XE g ) 2 } , 

(23) 

where the HOMO-LUMO exchange integral Kq = 
(<p _<t> 0+ \Vee\<Po+<Po-) C/[H]-(2i?)- 1 +C(exp.) reduces 
asymptotically to the atomic Hartree energy, and P is an 
(here) unspecified, yet smooth function of A. Note that 
the first two terms on the R.H.S. of Eq. 23 also follow 
from the single pole approximation to TDDFT excitation 
energies, and that the remainder describes corrections 
due to the coupling with higher KS excitations [84] . The 



corresponding eigenvector is = <^ +o-+C'(Ay / £Q. 
Hence we can split the asymptotic RPA response of H 2 
as follows: 

X A [H 2 ]=x A [H 2 ]+<5x A [H 2 ]+0(cxp.) , (24) 

where 

^[H 2 ](r,rV U ).-4^J^ , (25) 

and x A = (1 — X -^ w oo) _1 X° is the contribution from 
the other eigenvalues and eigenvectors of the matrix 
Mij^i(X). For large R, we can further write x A as the 
sum of the RPA responses of the H atoms, x A [Ha,b] ; and 
an inter-atomic correction [6] , Ax A [H 2 ] : 

X A [H 2 ]=x A [H A ]+X A [H B ] + Ax A [H 2 ] . (26) 

The atomic part 

X A [H]^(l-x°[H]A« cc )- 1 x°[H] (27) 

is formally equivalent to the RPA response of a spin- 
polarized H atom. 

Using the response functions as decomposed in Eqs. 
22, 24 and 26, the RPA correlation energy for stretched 
H 2 reads 



£ C RPA [H 2 ] = -/ 1 dA/ oo |HTrK{x A [H 2 ]M-x°[H 2 ](m)}] (28) 
= S C RPA [H A ]+S C RPA [H B ] + AS C RPA [H 2 ]+^ C RPA [H 2 ] . (29) 



Here £ RPA [H] is the RPA correlation energy of a 
spin-polarized H atom (associated with x A [H] — x°[H]), 
A£ RPA [H 2 ] comes from Ax A [H 2 ], and <5£ RPA [H 2 ] from 
<5x A [H 2 ] — (5x°[H 2 ]. Expanding Ax A [H 2 ] in Xv cc one can 
show, analogously to the case of interacting closed-shell 
systems [6], that the leading, second-order contribution 
to Ai? RPA [H 2 ] recovers the van der Waals interaction be- 



tween the H atoms, i.e. 

A£ RPA [H 2 ] ~ -C RPA [H]iT 6 , (30) 

with the atomic Ce coefficient obtained within the RPA. 
We next consider the HOMO-LUMO part of the re- 
sponse, writing <S£ C RPA [H 2 ] = /„* dX <5J7 RPA [H 2 ](A), where 



JC/ RPA [H 2 ](A) = - f °°f Tr[v ee {S X x M-6x (iu)}} 

= K ° (ofl) - !) (31) 
-^o({l + 4A|, + ^) +0 (A 2 )}- 1/2 -l) . 



r 



The integration over A then yields asymptotically: The term SE^ PA is associated with the static correlation 

<5£ RPA [H 2 ] = -K + O{^K^W g ) (32) 
~ -t/[H] + (2R)- 1 +C(exp.) . (33) 
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due to the (nearly) degenerate HOMO and LUMO KS 
states. Adding £ X [H 2 ] ~ 2£ X [H] + U[H] — (2R)- 1 + 
0(exp.) to £; RPA [H 2 ] we last get for the RPA XC energy 

£ RPA [H 2 ] ~ 2£ RPA [H] - Cf PA R- 6 . (34) 

The kinetic, electron-nucleus, nucleus- nucleus and 
Hartree components of the total energy of H 2 are eas- 
ily shown to also approach those of the free atoms for 
R — ► oo. Hence the RPA obeys the expected result 
lim^oo E PPA [H 2 ] = 2£ RPA [H] . ' 

Several remarks are in order: 
Leading R- dependence of E§ PA [H 2 ] for R — ► oo: Equa- 
tion 34 states that the C PPA R~ 6 van der Waals term is 
the leading correction to the asymtpotic RPA XC energy. 
This finding rests on the result of Eq. 33 that the static 
correlation term SE PPA [R 2 ] follows the (2R)^ 1 behav- 
ior of K (R), i.e. that it contains no multipolc terms of 
higher power than R- 1 up to i?~ 6 . The latter holds if (i) 
the HOMO and LUMO are represented by linear combi- 
nations of s-like atomic functions, as is appropriate for 
large R and as we had assumed. Then higher order terms 
in Kq(R) decay in fact exponentially with i?, as can be 
seen from a multipole decomposition of Kq. A further 
co nditio n is that (ii) the HOMO-LUMO gap (and thus 
\J KqE 9 in Eq. 32) decays exponentially, as is expected 
for Kohn-Sham states (as opposed to a Hartree-Fock cal- 
culation, where E g oc R^ 1 ) and which we have verified 
numerically in the range R = 4 ... 10 bohr. Of course, for 
the bond lengths considered here the van der Waals term 
is in fact marginal: for instance, C RPA i?~ 6 ~ 8 meV and 
~ 0.1 meV at R = 5 and 10 bohr, respectively, i.e. more 
than an order of magnitude smaller than the total RPA 
errors given in Tab. Ill (using C RPA ~ 4.6 a.u., calcu- 
lated from the atomic dipole polarizability) . Clearly, the 
RPA dissociation energy curve is still dominated in the 
range R = 4 ... 10 bohr by the static correlation term of 
Eqs. 32 and 33, as we further discuss in the next para- 
graph. 

Repulsion at intermediate R and role of self- 
consistency: From Eq. 31 we can interpret the ratio 
a RPA (A) = K E g /n PPA (X) < K as a correction to 
the exact asymptotic adiabatic connection (Fig. 2) that 
decays exponentially with R — > oo, turning on static 
correlation. a RPA (A) yields a positive 0{^/ KqE 9 ) con- 
tribution to the RPA exchange-correlation energy (see 
Eq. 33). While this contribution is expected to die out 
exponentially like E g , it may still be significant around 
R = 10 bohr, showing up as a spuriously repulsive dis- 
sociation curve. Our E g <~ 10~ 2 cV at R = 10 bohr 
is indeed compatible with the <~ 0.2 eV error we find 
from our RPA calculation at this bond length (this esti- 
mate follows from the single transition model discussed 
at the end of this appendix). We cannot rule out that in 
a self-consistent treatment E g decays (sufficiently) more 
rapidly compared to our present non-selfconsistcnt cal- 
culation. 

Behavior of E^ F X \H 2 \ for R — ► oo: For the exact 
exchange kernel [46] the same analysis of the molecular 



correlation energy as for the RPA can be carried through 
with A replaced by A/2. For R — > oo, the static correla- 
tion term SE PPA+X [H 2 ] C/[H] + (2.R)- 1 + C(cxp.) re- 
mains the same as in the RPA. However the atomic terms 
in Eq. 29 do not vanish, in contrast to the RPA+X corre- 
lation energy of spin-polarized free H atoms. Indeed, for 
a free H atom the spin-density responses Xa^ PA+X [H] 
and Xff CT '[H] arc identical, as is easily seen from the 
spin-resolved Dyson equation (see e.g. Ref. [77]). Thus 
£ RPA+X [H] = and £ RPA+X [H] = -0.5 a.u. However, in 
the spin-compensated stretched H 2 both the spin-up and 
the spin-down (nonintcracting) KS electrons are found 
with a 50% chance on either nucleus. Thus what en- 
ters as the atomic term x°[H] in Eqs. 22 and 27 cor- 
responds in fact to the KS spin-response of a fictitious 
spin-compensated H atom (denoted H') with half occu- 
pied Is t and Is J. states. The corresponding RPA+X 
correlation energy appearing on the R.H.S. of Eq. 29 is 
therefore £ RPA+X [H'] ~ £ RPA [H]/2 (estimated to second 
order in Aw C c)- Consequently, limij;__ +00 (i? t RPA+x [H 2 ] — 
£ RPA+X [H]) ~ £ RPA [H] < 0. In an actual calculation 
we indeed find that the RPA+X potential energy curve 
drops below zero beyond i? ~ 7 bohr (see Fig. 5). 

The RPA yields £ RPA [H'] = £ RPA [H], i.e. the distinc- 
tion between H' and H does not matter. This (spurious) 
RPA self-correlation energy (£ RPA [H] w -23 mHa per 
atom, using the exact tin) appears for the free atoms 
and the stretched H 2 and thus cancels out. In a spin- 
polarized one-electron system, the exchange kernel fully 
cancels the Coulomb kernel, leading to £? RPA+X [H] = 0. 
In a spin-compensated two-electron system like H 2 , corre- 
lation should arise only from the interaction of opposite- 
spin electrons. For the RPA+X kernel it is easy to see 
that there is zero like-spin correlation up to second or- 
der in Xv cc (although non-zero in higher orders). The 
fact that the dynamical correlation associated with \ x 
does not vanish for the RPA+X kernel in the infinitely 
stretched H 2 may then be seen as a failure to suppress 
opposite-spin correlation between the electrons on either 
atomic site. 

Had we treated the stretched H 2 in a spin-polarized 
KS scheme, the two electrons would be localized with 
opposite spin on either nucleus right from the beginning. 
While we have not performed this calculation, we ex- 
pect from our above discussion that both the RPA and 
RPA+X yield lim^^ £ to t[H 2 ] = 2E tot [H] in this case. 

Model based on the HOMO-LUMO transition only: If 
we include in the response only the HOMO and LUMO 
KS states, we get a "minimal" , two-state model of H 2 in 
which all effects of higher transitions are ignored. Then 
S C [H 2 ] in Eq. 29 is given by just the static correlation 
term SE C [R 2 }: 

E C [U 2 ] ~ %(^i + 2ga-i )-jr 

~ -K + v / 2K E g / K , fori?^oo , 

(35) 

for the RPA (k = 2) as well as the RPA+X (k = 1). 
The asymptotic independence of this minimal E C [H 2 ] is 
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again that given by Eq. 33. Hence within the minimal 
model also the RPA+X, like the RPA, correctly yields 
the total energy of the dissociated H2 as that of two free 
H atoms. Indeed, inspection of Eq. 31 shows that both 
RPA and RPA+X recover the exact adiabatic connection 
for R -> 00, i.e. that SU C (X > 0) = -J7[H]. Note however 
that due to the closing HOMO-LUMO gap for R -» 00 



the initial slope dU* PA+x (\)/d\\ x=0 = -K$/E g even- 
tually diverges, as does the the GL2 correlation energy. 
Independently of our work, an analogous minimal model 
has been recently obtained by van Leeuwen and cowork- 
ers, studying total energy functional based on Green 
functions [85]. 
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